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i.  Objectives: 

This  report  investigates  the  use  of  Bi-spectral  related 
techniques  to  extract  detection/classification  clues  from  time- 
frequency  representations. 

Earlier  results  have  indicated  that  1.5-D  spectral  techniques  (a 
degenerate  version  of  the  Bi-spectrum)  has  potential  Signal  to 
Noise  Ratio  (SNR)  gain  over  conventional  techniques.  This  is 
partially  due  to  the  rejection  of  Gaussian  like  perturbations  by 
the  cumulant  based  techniques.  The  third  order  moment  for  Gaussian 
zero  mean  random  processes  is  essentially  zero.  It  is  assumed  that 
i)  the  independence  of  IPS  from  analytic  signal  representations 
and  ii)  the  Gaussian  character  of  the  noise  causes  minimal 
distortion  of  the  output  variables,  hence  permit  signal 
characterization  at  low  to  moderate  SNR's. 

The  performance  is  demonstrated  on  i)  synthetic  signals  and  ii)  on 
segment (s)  of  the  NOSC  data  set. 


1.  Introduction: 

This  report  examines  the  instantaneous  power  spectrum  (IPS) 
and  the  cumulant  based  version  of  IPS  (CUM)  to  address  the  question 
whether  or  not  these  techniques  offer  possible  improvements 
relative  to  conventional  spectral  processing  (i.e.  spectrogram, 
Lofargram,  periodogram) . 

Section  2  provides  the  mathematical  definition  for  IPS,  CUM, 

and  LOFAR  while  section  3  allows  a  theoretical  comparison  of  these 

techniques.   The   comparison   is   performed   using   a   deflection 

criterion  (i.e.  maximum  signal  magnitude  versus  standard  deviation 

of  the  detector  output  under  noise  only  conditions) . 

Section  4  contains  simulation  results,  while  section  5 
provides  actual  SONAR  data  results. 

Section  6,  and  7  contain  the  conclusion  and  program  listings, 
respectively.   Section   8  provides  the  list  of  references. 


2.  IPS,  CUM  and  LOFAR 

2.1  IPS: 

IPS  is  a  member  of  the  Cohen's  class  [1],  and  is  also  known  as 
the  real  part  of  the  Rihaczec  distribution.  Our  particular 
implementation  differs  in  that  filtering  is  employed  as  the  data  is 
processed  (equation  1)  .  Earlier  results  have  shown  that  in  contrast 
to  the  Wigner-Ville  Distribution  (WVD)  no  spectral  cross  modulation 
terms  are  generated.  However,  it  has  wider  spectral  peaks  and 
superimposed  spectral  auto  modulation  terms  [2].  These  auto  terms 
do  not  interfere  with  the  interpretation  of  T-F  surfaces,  but  can 
help  in  the  detection  of  a  weak  stable  spectral  line  component.  The 
digital  implementation  of  IPS  is  given  by 

IPS(n,k)    =-Y?J-2N1/2  (x{n)x*{n-m)    +  x*  {n)  x(n+m)  )  w{m)   ex.p-j2xkin/N 

eq.  1 
where  k  =  frequency  index,   n  =  time  index,   N  =  the  length  of  data 
used,  m  =  shift  parameter  and  w(  )  =  lag  window. 
Some  earlier  results  are  published  in  [2,3,4]. 

2.2  CUM: 

In  general  the  Bi-spectrum  is  defined  as 

where  the  cumulant  is  given  by  [5] 

Cx(l±,l2)    =  E[X*  (n)  X{n+lJ  X(n+12)  )  . 


A  degenerate  version  of  the  cumulant  is  given  by 


CX(0,12)    =  E(X*(n)X(n)X(n  +  l2)  )    =  E(\X(n)  \2  X(n  +  12)  ) 


When  replacing  the  expected  value  with  the  instantaneous  value 
the  degenerate  Bi-spectrum,  also  called  the  \\   D  spectrum,  becomes 

The  cumulant  version  of  IPS  is  given  by 

CUM{n,k)    ^ly^^2'1    {\x{n)  \2x'  (n-m)    +  \x{n)  \2x{n+m)  )  w{m)    ex.p-j2nkm/N 

eq.  2 
Cross  terms  at  locations  where  the  true  spectrum  has  no  support 
are  not  created.  Also  contrary  to  the  conventional  Bi-spectrum  none 
of  the  spectrally  related  components  are  suppressed. 
In  part  of  the  work  in  [6],  the  approach  of  Petropulu  [7]  which 
uses  an  instantaneous  higher  order  moment  slice  was  utilized.  This 
approach  works  well,  but  does  not  permit  separation  of  spectral 
components  having  the  same  spectral  dynamics.  That  is,  all 
stationary  spectral  components  (i.e.  stable  line  components)  map  to 
the  same  detection  cell,  making  this  technique  useless  in 
scenarios  with  several  similarly  behaving  spectral  components. 

2.3  LOFAR: 

The  LOFAR  algorithm  uses  a  Hamming  windowed  Fast  Fourier 
Transform  (FFT)  and  displays  the  magnitude  of  the  transform. 


To  keep  T-F  dimensions  relative  small,  no  zero  padding  is  used 

Px(n,k)    =  |  r^xl^e' 

eq.  3 
The  IPS,  CUM,  and  LOFAR  surfaces  are  averaged  (in  magnitude)  over 
some  segments  of  the  time  axis  leaving  only  a  display  of  spectral 
magnitude  versus  frequency.  For  each  segment  one  spectral  line  is 
created  resulting  again  in  a  time  frequency  display.  The  averaging 
technigue  is  also  known  as  power  or  incoherent  averaging. 

3.  Theoretical  results  for  averaged  surfaces. 

The  analysis  in  this  section  is  done  by  disregarding  any 
window  effects  and  assuming  that  under  the  noise  only  hypothesis 
(HQ)  Gaussian  white  noise  and  under  the  signal  only  hypothesis  (H^ 
a  deterministic  signal  of  the  form  x(n)=A  cos  (27rkQn/N)  is  present. 
The  value  of  kQ,  n,  and  N  are  assumed  to  be  integers  which  ensures 
that  the  spectrum  peaks  at  a  bin  center  of  an  N-point  Fourier 
transform.  In  all  detection  schemes  (averaged  IPS,  averaged  CUM, 
and  averaged  LOFAR)  magnitudes  are  averaged  since  they  are 
numerically  more  conveniently  computed  then  the  magnitude  squared 
values. 

The  comparison  is  done  for  normalized  versions  of 
I(k)  =  Sn  |lPS(n,k)  |  , 
C(k)  =  Sn  |cUM(n,k)  |  ,   and 

P(k)  =  Zn  |windowed  Fourier  transform  {x(n1,k)}|  ;  where  n1  is  the 
time  counted  from  reference  point  n. 
In  what  follows  the  time  (n)  and  frequency  (k)  dependency  is 


usually  suppressed 


3.1  Averaged  IPS. 


x(n) > 


IPS(n,k) 


real 


--J 


i  i 


>  1/M  Z 


y  z  w 

where  x(n)  ~  N(0,a2)  i.i.d.  and  M   is  the  appropriate  number  of 
degrees  of  freedom. 


Under  H0  : 

Under  noise  only  condition  (HQ)  the  random  variable  Y  denoted  by  YQ 

has  the  following  moments  [6]: 

E{Y0)=  a2    ,  E{Y0)2  =  [(N+5)/2]  a4  ,  var(YQ)  =  [(N+3)/2]  a4  . 

For  the  transformation   Z  =  SQRT(Y2)  , 

the  resulting  probability  density  function  is  given  by 

f2(z)  =  (fY(z)  +  fY(-z))  U(Z)  . 

Making  the  assumption  that  for  large  enough  number  of  data  points 

the  random  variable  Y  is  approximately  Gaussian  distributed  we  can 

derive  the  moments  for  Z  as 

E(Z0)  =  SQRT{  (N+3)/167r)  a2  e"(2/(N+3)) 

E(Z02)  =   E{Y0)2  =  ((N+5)/2)  a4 

var(Z0)  =  [(N+5)/2  -({(N+3)  e"(/;/(N+3l>  }/16tt)  ]  a4 

var(W0)  ss  l/M   var(Z0)   (for  large  N) 

w  (1/M) 0.48  N  a4. 
A  typical  surface  has  about  N/3  degrees  of  freedom  [6],  providing 
a  standard  deviation  (under  the  H0  hypothesis)  of 


<7(W0) 


1.2    O' 


Under  H,      (at   spectral    location   ±   kQ) 
Z1    =    |  Y1  |     =    (NA2)     cos2(27rk0n/N) 
W1    =    (NA2)/M      En   cos2(27rk0n/N) 
»    (NA2)/2. 


The   deflection   DpsI    for   IPS   becomes 


2\     _ 


.2  /„2> 


Dp„    =   W1/a(Wn)    =    (NAV2)/(1.2    o<)    =    (N/2.4)  (AV^) 


3.2   Averaged  CUM. 


x(n) > 

CUM 

C 

real 

i      i 
I      i 

1/M   Zn 

— > 

Under  H0: 

Y 

Z 

W 

E(Cn) 


o  +  jo 


r2\   3 


Z0     " 


E(C0)2   -  var(C0) 

=(6N+54)/4(a')°        var(YQ)     ;     [6] 
(since   real   part   of   cum   is    zero) . 

Yo    I 
fz(z)=    2/(27raY2)1/2   exp- (z2/ (2ay2)  )    U(z) 

E(Z0)    =    SQRT(2/tt)     aY 

E(Z0)2   =    3/2    aY2 

var(Z0)    =    (3/2    -    2/tt)    ay2 

E(W0)     =    E(Z0)     =    SQRT(2/tt)     ay 

E(W„)2    =    1/M    (E(ZQ)2   +(M-l)E2(Zn)  } 


var(W0)  =  l/M  {  (  3/2-2/w)  (6N+54  )  /4  }  a6 

a(W0)  =  SQRT{1/M  (3/2-2/7T)  (6N+54)  }/2  a3. 

A  typical  surface  has  about  N/3   degrees  of  freedom  [6]  ,  so  the 

standard  deviation  becomes 

a(Wn)  «  SQRT{3/N  (3/2-2/7T)  (6N+54)  }/2  a3. 


Under  H,  (at  spectral   location  ±  k0  )  : 
x(n)  =  A  cos  (27rk0n/N) 

Y1  =  NA3/2  cos3(27rkn/N)      @  k  =  ±  kQ 
Z1  =  NA3/2  |  cos3(27rkn/N)  |     @  k  =  ±  kQ 
From  numerical  evaluations  we  know  that 
bounded  by 


the  summation  is  lower 


-y^'1  I  (cos3  (2-Kln/M)  )  I  >  0  .  424 


Hence    Z,    >    0.424    NA3/2 


The   deflection   DCUM   becomes 


CUM 


(0.424    NA3/2     )SQRT(M)/{SQRT[  (3/2-2/7T)  (6N+54)/4]     a3} 
«    0.1075    N     (A/a)3    . 


3.3    LOFAR: 


x(n) ; 


X(k) 


i     i 
i     i 


l/M    Sn  > 

W 


Under   H0: 


X(k)=    Xr(k)     +    j     X,(k) 


=   2  _  rv+N"1     x(m)    exp-j  (27rkm/N)  ;    where   n    is   the   reference   time 
E(Xr(k))=    E(Xr)=    E(Xi(k))=    E(Xi)=    0 
E(Xr)2   =   var(Xr)    =   var(X.)    -   a2   N/2 
Z    =    SQRT(Xr2   +    X-2    ) 
E(Z0)     =SQRT(7T/2)     SQRT(N/2)     O 
var(Z0)    =    (2-7r/2)a2   N/2 
E(W„)    =    E(Z„)     =SQRT(7T/2)     SQRT(N/2)     a 
E(W0)2    =1/M    (E(Z0)2   +     (M-1)E2(Z0)} 
var(W0)    =    1/M  var(ZQ)    =    (1/M)     (2-7T/2)  a2  (N/2) 
a(W0)    =    SQRT{(1/M)     (2-W/2)     (N/2)}    a 
Under  H.,    (at   spectral      location   ±   k0  )  : 
XJk)    =A  N/2        @   k  =  ±   k0 
Z1    =  A  N/2  @    k  =   ±    k0 

E(W1)     =    E(Z,)     =   A   N/2 
The   deflection   D. nc.D   becomes 

LOFAR 

Dlofar    =    (AN/2)    /{(1/M)  (2-7r/2)a2(N/2)}1/2 
=   N1/2    [3    A2   /(  (4-7T)a2)  ]1/2    . 

3.4  Deflection  ratio  comparison: 

We  can  now  plot  the  deflection  ratio,  that  is  compare  and 
determine   which   detector   is   more   advantageous   under   which 
conditions. 
Comparing  D,ps  and  DL0FAR  : 

d.ps   >   DLOFAR 
simplifies  to 

N  >  20.13  (a2/A2) 

Comparing   DCUM    and    DL0FAR    : 


10 


CUM    >    DLOFAR 


simplifies  to 

N  >  302.4  (aVA4)  . 


We  notice  that  the  processing  length  variable  N,  for  IPS  relative 
to  LOFAR  has  a  (SNR)"1  and  for  CUM  relative  to  LOFAR  has  a  (SNR)"2 
dependency. 


4.  Simulation  results  for  sinusoid  in  white  Gaussian  noise. 

Ten  test  sets,  of  length  4096,  consisting  of  sinusoidal 
signals  embedded  in  white  Gaussian  noise  are  processed  by  IPS,  CUM 
and  LOFAR.  Each  test  set  contains  seven  (7)  sinusoids  whose  SNR 
goes  from  -3  dB  to  -21  dB  in  3  dB  steps.  The  decrease  is  monotonic 
relative  to  the  spectral  locations  so  that  the  SNR  decreases  in  3 
dB  steps  as  the  frequency  increases.  The  sinusoidal  frequencies 
are  fixed  in  all  test  sets  at  digital  frequencies: 
(a)  }  =  {IIV3,  17%,  29V3,  37V3,  59V3,  77V3,  lll%}*27r/256  .  The  phases 
within  each  set  and  over  all  set  are  independently  distributed 
(i.e.  no  harmonic  relationship  between  different  spectral 
components  and  different  realizations) .  We  note  that  none  of  the 
sinusoids  has  an  integer  number  of  cycles  over  4096  data  points 
ensuring  all  signals  will  have  side  lobes  if  processed  by  an  FFT  of 
size  4096  or  less.  To  demonstrate  potential  gain  the  deflection 
ratio  as  defined  in  section  3  is  computed  for  each  sinusoid  and 
averaged  over  the  10  realizations. 

Figures  1,  2,  3,  4,  5  and  6  are  plots  of  the  outcomes  of  the 
simulation  for  128,  256,  512,  1024,  2048  and  4096  data  points, 


11 


respectively.  The  top  part  of  each  figure  superimposes  all  of  the 
10  realizations  while  the  bottom  part  shows  the  average  of  the  10 
realizations.  We  note  for  example,  that  the  weakest  sinusoid  (bin 
223-224  of  512  point  data  length)  clearly  shows  up  in  the  CUM  and 
IPS  representation  (figure  3. A)  but  is  not  detectable  in  the  LOFAR 
representation  (figure  3.B).  This  illustrates  that  the  alternative 
representations  (i.e.  averaged  IPS  or  averaged  CUM)  can  complement 
the  traditional  (averaged  LOFAR)  representation.  Figures  7-12 
provide  plots  of  the  experimental  deflection  ratios  versus  SNR  in 
the  top  half  of  the  figures.  The  solid,  dashed,  and  dotted  lines 
represent  the  experimental  deflection  for  CUM,  IPS,  and  LOFAR, 
respectively.  The  bottom  halves  are  plots  of  the  difference  of 
deflection  ratios  (D,ps  -  DCUM  )  .  Figure  13  attempts  to  reconciliate 
theoretical  bounds  from  section  3  with  the  experimental  evidence 
from  this  section.  The  straight  solid  lines  demarcate  the 
theoretical  transitions  points  at  which  CUM  works  better  than  LOFAR 
(upper  straight  line)  and  where  IPS  works  better  than  LOFAR  (lower 
straight  line) .  We  mention  again  that  no  window  action  was 
accounted  for  in  the  derivation  of  these  lines.  Superimposed  are 
(zigzag  lines)  the  results  taken  from  figures  7  trough  12.  The 
upper  zigzag  line  shows  the  demarcation  line  where  the  experimental 
evidence  indicates  IPS  will  do  better  than  CUM  (in  terms  of 
deflection  ratios) .  The  lower  zigzag  line  consists  of  the  points, 
as  a  function  of  SNR  and  processing  length,  where  the  experimental 
performance  of  all  three  technigues  (IPS,  CUM  and  LOFAR)  tends  to 
coincide.  Above  this  line  the  CUM  and  IPS  technique  outperform  the 
LOFAR  approach  with  improvements  increasing  as  the  input  SNR  increases. 


12 


ipsav(x1 28,1,1 28,4) 


curnav(x128,1, 128,4) 


700 


600- 


average  of  above 


average  of  above 

i  iinn 

1  '-tUU 

1200 

., 

i           i           i 

_ 

1000 

1 

800 

- 

- 

600 

1 1 

400 

?nn 

\J 

i 

20 


40 


60 


Fig.    l.A 


L3 


lofarm(x1 28,1 28,1 281 28) 


35 
30 

average  of  above 

■              i              i 

25 

j 

20 

' 

15 

i 

10 

- 

- 

5 
n 

A/J 

U 

— i 1 i 

20 


40 


60 


Fig.    i.b 


14 


ipsav(x256, 1,256,4) 


curnav(x256,1,256,4) 


3500 


7000 


20  40  60  80  100  120 


20  40  60  80  100  120 


3000 


2500 


2000 


1500 


1000 


500 


average  of  above 


t 1 r 


uMW 


o 


i i 


20  40  60  80  100  120 


5000 

4500- 

4000- 

3500 

3000 

2500 

2000 

1500 


1000 


500 


average  of  above 


-i 1 r 


J 


w 


-I 1- 


i _ i_ 


20  40  60  80  100  120 


Fig.  2. A 


15 


lofarm(x256,256,256,256) 


20  40  60  80  100  120 


average  of  above 


20  40  60  80  100  120 


Fig.  2.B 


16 


ipsav(x51 2,1,51 2,4) 


14000 


12000- 


x  10ftumav(x51 2,1,51 2,4) 


50   100  150  200  250 


50  100  150  200  250 


12000 
10000 
8000 
6000 
4000 
2000 
0 


average  of  above 


^ 


t 1 r 


u 


Y^ 


W 


ty&k 


'iWW/yw 


-i i i 


X10   ai 

i/erage 

of  above 

1.8 

i       i       i       i       i 

1.6 

- 

1.4 

- 

1.2 

- 

1 

- 

0.8 

- 

0.6 

- 

0.4 

- 

- 

0.2 

K 

i  J 

mWmJ 

iWwW^ 

50   100  150  200  250 


50  100  150  200  250 


Fig.  3. A 


17 


lofarm(x5 1 2,5 1 2,5 1 2,5 1 2) 


160 
140 
120 
100 

80 

60- 

40- 

20 


-i r 


-i 1 r 


0rrrrfffT 


IJ  lLjiLJiIi 111 


l^yulL'k 


50      100     150    200    250 


140 

120 

100 
80- 
60- 
40- 
20- 


J*i 


average  of  above 


t r 


t r 


VUw 


U^ 


o 


-1 L 


50      100     150    200    250 


Fig.    3.B 


18 


x1(jfisav(x1 024,1,1024,4) 


4.5 

4 
3.5 

3 
2.5 

2- 
1.5- 

1 
0.5 


t r 


JJI 


iJil  L 


Imy.l  .njLjiJ«li 


J.L.^a^^j.^..  J-».iLli1  ■*■  A   ^U^JJI 


x1gftmav(x1 024, 1,1 024,4) 


100     200    300     400     500 


100    200     300    400     500 


4.5 


x  -|Q4  average  of  above 


4 
3.5  h 


•3  h 


2.5 


2 
1.5 

1 
0.5 

0 


1 r 


A* 


I 


WK"wrt»*Ww 


JlfWv^^Ar^v 


100     200    300     400     500 


x  -|Q4  average  of  above 


1 


0 


M-Ww 


W»^v^V^rvU^A^^^^^ 


Fig.    4. A 


100    200     300    400     500 

19 


lofarm(x1 024,1024,1 024,1 024) 


300 
250 
200 
150  - 

100- 
50 
0 


t 1 ( 1 r 


iL  i 


li^L.^i 


M 


Ijllu  i^lM^i^M- 


ltii.ji.Jk..   »ij.llt,L— ^.lLll.  . 


100  200  300  400  500 


300 


250- 
200- 
150- 
100- 
50 
0 


average  of  above 


-i 1 1 r 


pAYWW 


^mms^i^^ 


100    200    300    400    500 


*ig.    4.B 


20 


x1(jfeav(x20348, 1,2048,4) 


x  1qj5mav(x20348, 1,2048,4) 


18f 


3 

2.5 

2 

1.5 

1 

- 

0.5 

.JulJl 

200    400     600     800    1000 


0 


jiki  tiiiu  -^.f^+m. 


200  400  600  800  1000 


x  1Q4  average  of  above 


200  400  600  800  1000 

Pig.    5. A 


3.5 

3 
2.5 

2 
1.5 

1 
0.5 

0 


x  -|q5  average  of  above 


-| r 


I L 


m,^.X^iM»pnWI^iViJ»>A*^ 


200    400     600     800    1000 


21 


Iofarm(x2048,2048,2048,2048) 


600 
500 
400 
300 
200 
100 
0 


200  400  600  800  1000 


average  of  above 


600 
500 
400 
300 
200 
100 
0 


t r 


■WWirtWhft  hfH  V^Ui^frylidiM. 


J l_ 


200    400    600    800   1000 

Fig.    5.B 


22 


x  -|j(Ssav(x4096J, 4096,4) 


8 
71- 

6 
5 
4 
3 
2 
1 
0 


x  1ig§mav(x4096J,4096J4) 


500      1000     1500     2000 


500      1000     1500     2000 


x  ig5  average  of  above 


6- 
5- 

4- 


0 


m*hM*Jmmil***htm*M*to*mi**** 


14 

12 

10 

8 


x  igb  average  of  above 


500   1000  1500  2000 

Fig.    6. A 


0 


w+*tm*m*mi* m 


km*  *Hmmm*ti¥ 


1 


500      1000     1500     2000 


23 


Iofarm(x4096,4096,4096) 


1200 


1000 


800 


600 


400- 


200- 


I.  il    I    J  L    I    .1.  l.i  ill  J  Id  I  lij.j    jIlLkLl 


0 


t-  J~  I  ImJ-J...  —i..  I  -^  h.  i-  ..Ik.-  .X.  .    uLl    xJl,_  .Llft.lL&.k 


500      1000     1500     2000 


1200 


1000h 
800 
600 
400 
200 
0 


average  of  above 


ii 


I  L. 


i  .IL 


Ul  lij..    .ikiLLl 


J,.  . 

—  I  LJ.J ...  ^_— .k   ^  W.   1..  ,lk.-    .X.    .     Utl      <_»L..LiA.(L»..kl 


500      1000     1500     2000 

Fig.    6.B 


24 


data  size  128,  proc  length  128 

9 

i                 i                 i                 i                 i 

8 

-  \* 

7 

NX 

6 

\ 

•Ss 

o 

CD 

u4 
Q 

\ 

3 

'"•    ^^ 

2 

'"•■•■.     ^v, 

1 

'\          - 

n 

i                          i                         i                         i                          i 

3  4  5 

SNR  =  -3dB  times  index 


6 


7 


0.25 


data  size  256,  proc  length  256 


3 
SNR 


4  5 

-3dB  times  index 


7 


Dp21-Dq21 

Fig.    8 


_0 


o 

CD 

"5= 
CD 


data  size  512,  proc  length  512 

20 

NX 

■                i                i                i 

15 

- 

10 

- 

N. 

5 

i 

i                 i                 i                 i        '•■•■•. 

3 


4  5 

SNR  =  -3dB  times  index 


7 


0.1  F 


data  size  1024,  proc  length  1024 


3  4  5 

SNR  =  -3dB  times  index 


Dpi  1-Dql 1 

Fig.    10 


28 
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5.  Processing  Results  using  NOSC  Data  Base. 

In  this  section  the  data  set  supplied  by  NOSC  is  processed. 
Figure  14  is  a  LOFARGRAM  of  the  data  set  r2A_14  which  is  4,128,769 
data  points  long.  The  time  origin  for  all  grams  is  chosen  to  be  at 
the  top  of  each  figure.  Figure  15.  A  and  15.  B  are  LOFARGRAMS 
produced  at  NPS  to  verify  the  fidelity  of  the  data  transfer  from 
NOSC  to  NPS.  To  see  representative  spectral  plots  we  selected  the 
most  interesting  segment,  that  is  from  time  sample  1,024,001  to 
1,228,800.  This  corresponds  to  the  region  508  to  408  along  the  time 
(vertical)  axis  of  figure  15. A.  Initial  processing  using  IPS  and 
especially  CUM  showed  heavy  modulation  of  the  time  freguency 
displays.  This  is  caused  by  the  low  frequency  spikes  observable  on 
figure  14  and  15. A  (i.e.  strong  DC  like  pulsations  (about  5800 
sample  points  apart) .  In  the  top  part  of  figure  16  which  displays 
the  time  series,  used  in  the  experiments,  these  modulation  spikes 
are  clearly  seen.  To  minimize  the  unwanted  amplitude  modulation 
effects  on  the  time  frequency  plots,  the  data  is  soft  limited  by 
scaling  all  values  larger  than  +90  by  dividing  them  by  4  and  by 
scaling  all  values  less  than  -90  by  dividing  them  by  2.  The 
resulting  time  series  is  plotted  in  the  bottom  part  of  figure  16. 

The  contour  plots  of  the  remainder  of  this  section  show  now 
little  amplitude  modulation  allowing  the  approximation  to  gray 
scale  via  contour  plots  to  work.  On  a  color  monitor,  with  a 
sufficient  number  of  color  levels  of  quantization,  the  modulation 
may  not  be  bothersome.  Figure  17  displays  49  traces  (along  the  y- 
axis)  which  corresponds  to  49  individual  estimates  of  the  spectral 
density.  Closer  examination  reveals  some  details  on  the  IPS  plot 
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not  so  easily  seen  on  the  lower  LOFARGRAM  plot.  In  particular  we 
refer  to  3  areas: 

1)  tick  mark  38-48  (vertical   axis)  about  bin  700  (horizontal 
axis) , 

2)  tick  mark   1-33  (vertical   axis)  about  bin  850  (horizontal 
axis) ,   and 

3)  tick  mark   1-20  (vertical   axis)  about  bin  1300  (horizontal 
axis) . 

A  word  of  caution  however  when  using  the  simple  minded  contour 
(MATLAB)  plots:  peaks  of  identical  height  but  of  different  width 
will  appear  differently.  In  appendix  A  the  effects  of  choosing 
different  contour  levels  is  illustrated.  However  for  all  plots  the 
best  choice  of  level  was  manually  selected  by  examining  several 
plots  differing  in  the  number  of  levels  selected  and  then  retaining 
the  one  with  the  best  features. 

Figure  18. A  and  18. B  are  displays  of  IPS,  CUM  and  LOFAR  with 
the  lines  created  512  time  data  points  apart.  Figure  18. A  (as  well 
as  19. A  and  20. A)  displays  spectral  location  51  through  1500  while 
figure  18.  B  (as  well  as  19. B  and  20.  B)  displays  spectral  location 
351  through  1500.  This  allows  better  interpretation  of  the 
different  spectral  regions. 

Figures  18,  19  and  20  have  the  same  processing  parameters, 
they  differ  in  the  number  of  spectral  lines  computed.  Figure  18 
uses  a  shift  of  512  (i.e.  8:1  overlap)  resulting  in  390  lines, 
figure  19  uses  a  shift  of  1024  (i.e.  4:1  overlap)  resulting  in  197 
lines,  while  figure  20  uses  a  shift  of  4096  (i.e.  no  overlap) 
resulting  in  49  spectral  lines. 
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Figure  21  illustrates  performance  as  the  processing  length  is 
increased  to  8102  (shift  8192,  no  overlap)  providing  25  spectral 
lines  (bin  51  through  3300  are  displayed) . 

6.  Conclusion  and  recommendation. 

As  illustrated  in  section  3  and  4  processing  gain  can  be 
obtained  for  IPS  and/or  CUM  relative  to  the  LOFAR  processing 
technique.  This  is  a  function  of  signal  and  processing  parameters. 
It  may  allow  glimpses  of  targets  which  otherwise  maybe  hidden. 
However,  the  processing  gain  when  achievable  is  bought  at  the 
increase  in  processing  cost.  IPS  and  CUM  are  actually  tools 
designed  more  for  transient  type  signals  but  as  shown  can  be  used 
to  obtain  gain  on  stable  narrow  band  signals.  Simplified 
processing  was  not  attempted  but  deserve  future  examination.  The 
analysis  in  section  3  is  done  totally  disregarding  lag  and  data 
windows  which  typically  are  used.  Future  analysis  should  address 
the  influence  of  the  window  function  on  processing.  Extensions  to 
extract  phase  information  is  not  attempted  but  should  be  pursued. 
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ipsav(rf,1, 4096,32)  bin  351:1500,  shft  1024 


cumav(rf,1,4096,32)  bin  351:1500,  shft  1024 


150 

100 

50 


7^ 


a. 


i    . 


i..;*.:w;kf 


F I'M'    1'|,L     "^         >'*r 


n  .     .     ■  ft. 

*    .»L  > 


.':*l:fi'  •?• 


i 

•r. 


ft; 


■   " 


■  ■•  a 


r-v '-,' 


.(i*.  • 


^        'I  '  •  1  r .'  1   'I  *  ■  '■      J  tf     ».'r      *  |J  1  i  I^M  I'-J     I  .  <T.%i  '"A     •  H  .1   -#4  .  I 


»■'   ' J i_a_ t ,:;    ''  ,   VI ■    '•-■■'    *■■ i    <■  ~-l_^. >V^     ■  Ji  n.itts. ,___ 


200 


400 


600 


800 


1000 


lofarrnfrf, 4096,4096,4096)  bin  351:1500,  shft  1024 


■■Itfl  :U    -    >•.   v    'i  mH 


150 

100 

50 


II 


.    4      » 


IF*  *f  >  Jf  ?.<  <  "d:  ;'£  IT  ill  fiV 


■v-i. 


1" .     1", 


200  400  600  800  1000 

soft  clipped  data,  gray  levels  top  to  bottom  4,4,4 


Fig.    19. B 


43 


ipsav(rf,  1,4096,32)  bin  51:1500,  snft  4096 


t — n — I n — r~~^ — '<r<  "  •-■ — ~ ' — m* — •-*-■ — r  prim — fr: — i. -.ii-  ; 


40 
30 
20 
10 


>, 


-i 


•  ivf 


Iff 


■^ 


»  >  i- 1 


«♦  ,/m  , .;  ..•!.,  >;h  4*  III"  Vs  1  iu..,u\*V*v"J 


200 


400 


600 


800 


1000 


1200 


1400 


cumav(rf,  1,4096,32)  bin  51:1500,  shft  4096 


40 
30 
20 
10 


Mil   J  i  fi 

t"'»"ji  Ii"' 


!-■► 

♦ 


"*'|Hji*    »•      m'H  -►!-«-•'  V^fc  *|    -l^iiltM'    PVWf*    WtV 


200 


400 


600 


800 


1000 


1200 


1400 


lofarm(rf,4096,4096,4096)  bin  51:1500,  shft  4096 


40 
30 
20 
10 


•A 

I 


■•MM 


200  400  600  800  1000         1200 

soft  clipped  data,  gray  levels  top  to  bottom  6,7,6 


1400 


Fig.    20. A 


44 


ipsav(lt  1,4096,32)  bin  351:1500,  shft  4096 
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7.  Program  Listings: 

All  algorithms  are  implemented  (including  their  hard  copy 
functions)  using  MATLAB  software  (The  Mathworks  Inc.) .  The  programs 
can  be  executed  on  MS-DOS  personal  computers  (i.e.  386  +  co- 
processor or  486  based  IBM  compatible  PC's)  or  on  SUN  work 
stations.  Due  to  the  long  data  sequence  and  large  T-F  surfaces  most 
of  the  work  is  done  on  SUN  work  stations.  The  hard  copies  provided 
are  obtained  by  using  the  contour  plotting  function  of  MATLAB.  One 
note  of  caution  is  appropriate  when  interpreting  the  T-F  plots.  Due 
to  the  way  contour  lines  are  drawn,  equal  height  spectral  peaks  can 
appear  different  if  they  differ  in  T-F  spread,  that  is  a  narrow 
peak  will  be  darker  than  a  wide  peak  even  though  these  peaks  are  of 
the  same  height.  This  is  due  to  the  density  of  contour  lines  which 
is  higher  for  narrow  peaks  than  it  is  for  wide  peaks.  Better 
results  are  obtained  when  displaying  the  data  on  gray  scale  or 
color  scale  output  devices.  The  new  4.0  MATLAB  version  supports  the 
function  'image'  which  on  a  gray  scale  or  color  monitor  will 
improve  the  readability  of  LOFARGRAMs. 

The  color  output  (using  MATLAB  3.5i)  on  a  SUN  SPARC-2  work 
station  shows  superior  detail  when  compared  to  the  hard  copies 
supplied  in  this  report.  The  hard  copies  are  obtained  using  a 
screen  dump  routine  which  preserves  less  detail  than  a  postscript 
derived  output.  The  postscript  output  is  not  used  due  to  the  length 
of  time  involved  generating  it  on  the  processing  system. 
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>>  atypical   setup  for  obtaining  IPS  and  CUM  like  surfaces 

>>  %for  i=start:finish 

>>  *pCi,:)HpsavCrf((i-1)*shft+1:Ci-1)*shft+N),1,N,step); 

>>  *qCi,:)=cumavCrf(Ci-1)*shft+1:Ci-1)*shft+N),1,N,step); 

>>  %   N=  proc  length,  step=stepsize,  N/step=number  of  terms  averaged  to 

>>  ^create  one  spectral  line,  shft  =  amount  of  data  shifted  to  create 

>>  %the  next  surface  which  will  be  averaged  to  create  one  spectral 

>>  %line  at  location  i 

>>  %end 

>> 

>>  atypical    setup  for  obtaining   LOFAR  like  surfaces 

>>  %for  i=start:finish 

>>  %po(i,:)=lofarm(rf((i-1)*shft+1:(i-1)*shft+N),id,st,ti;)j 

>>  %id=number  of  data  points  in  transform,  st=stepsize, 

>>  %tl=transform  length 

>>  %  Note:    rf  =  input  sequence, 

>>  %  Note:   in  all    run  executed  id=st=tl 

>>  %  So  the  overlap  was  controlled  with  shft 

>> 

>> 
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%function  [P, freqindex] =ipsav (data, wintype, winlen, step) ; 

%This  function  will  calculate  a  single  averaged  Instantan.  Power  Spectral  (IPS) 

%  line. 

%The  IPS  surface  characteristics  are  determined  by  the  selection  of  window 

%type  (wintype) ,  window  length  (winlen)  and  the  distance  that  the  window 

%is  moved  through  the  data  sequence  (step) . 

s- 

o 

%The  P  matrix  plots  only  positive  frequencies  (in  magnitude) .   The 

%outputs  timeindex  and  freqindex  can  be  used  in  plots  to  interpret  the 

%results  . 

%The  inputs  are: 

%data     -  The  input  data  string  row  vector 

Iwintype:   '0'  Rectangular  Window 

%  '  1'  Hamming  Window 

%winlen   -  The  desired  width  of  the  window,  normally  half  of  the  siglen 

%step     -  Time  step  desired,  normally  ' 1'  or  a  multiple  of  '2' 

%See  also  IPSSURF,  IPSLOFAR 

function  [P, freqindex] =ipsav (data, wintype, winlen,  step) 
[datarows, datacolumns] =size (data) ; 
if  datarows  ~=1 
data=data . ' ; 
end 

siglen=length (data)  ; 
if  wmtype==0 

wm=ones  (wmlen-1,  1)  ; 

elseif  wintype==l 

win=hamming (winlen-1) ; 
end 

W= [win (winlen/2 :-l : 1) ]  ; 

x= [zeros (1, winlen)  data  zeros (1, winlen) ] ; 
P= zeros (1, winlen/2) ; 

for  n=winlen+l : step : siglen+winlen-step+1 

Xm= [con j  (x (n:-l :n-  (winlen/ 2-1) ) )  . '  x (n:n+ (winlen/2-1) ) .'  ]  ; 

Xn=[x (n) ; con] (x (n) ) ]  ; 

product= ( (Xm*Xn)  . *W)  . '  ; 

product= [product  0  con j (product (winlen/2 :-l : 2) )]  ; 

ptemp=f ftshift (real ( . 5*f ft (product)  )  )  ; 
P=P+abs (ptemp (winlen/2  +  1 : winlen)  )  ; 
end 

[prow,pcolumn] =size (P) ; 
freqindex= [ 0 :pcolumn-l ] ; 
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%function  [P,  f  reqmdex,  timeindex]  =cumav  (data,  wintype,  winlen,  step)  ; 

%This  function  will  calculate  the  1.5  D  Spectral  surface. 

%this  surface  consists  of  averaged  (shorter)  cumulant  based  surfaces 

%in  a  power  or  magnitude  averaged  sense 

%The  1.5  D  surface  characteristics  are  determined  by  the  selection  of  window 

%type  (wintype) ,  window  length  (winlen)  and  the  distance  that  the  window 

%is  moved  through  the  data  sequence  (step) . 

%The   surface  is  placed  sequentially  in  the  P  matrix  for  display. 

%The  P  matrix  plots  only  the  positive  half  of  the  spectral  plane.   The 

%outputs  timeindex  and  freqindex  can  be  used  in  plots  to  interpret  the 

%results . 

%The  inputs  are: 

%data     -  The  input  data  string  must  be  in  row  vector  form 

%wintype:   '0'  Rectangular  Window 

%  ' 1'  Hamming  Window 

%winlen   -  The  desired  width  of  the  window,  normally  half  of  the  siglen 

%step     -  Time  step  desired,  normally  ' 1'  or  a  multiple  of  '2' 

%See  also  ONESURF,  ONELOFAR 


function  [P,  freqindex] =cumav (data, wintype, winlen, step) 

[datarows, datacolumns] =size (data) ; 

if  datarows~=l 

data=data . ' ; 

end 

siglen=length (data) ; 
if  wintype==0 

win=ones (winlen-1, 1) ; 
elseif  wintype==l 
wm=hamming  (winlen-1)  ; 
end 

W=[win  (wmlen/2  :-l:l)  ]  ; 

x= [zeros (1,  winlen)    data   zeros (1, winlen) ] ; 
P=zeros (l,winlen/2) ; 
for  n=winlen+l : step : siglen+winlen-step+1 

Xn=[abs  (x(n)  )  "2    ;  abs  (x  (n)  )  "2  ]  ; 

Xm= [con] (x (n:-l : n-  ( winlen /2-1) ) )  . '  x (n:n+ ( winlen /2-1)  ) .'  ] ; 

product= ( (Xm*Xn) . *W) . ' ; 

product= [product  0  con j (product (winlen/2 :-l :2)  )]  ; 

ptemp=f ftshift (real ( . 5*f ft (product) ) )  ; 
P=P+abs  (ptemp  (winlen/2  +  1  :wmlen)  )  ; 

end 

[prow, pcolumn] =size  (P) ; 


f reqindex= [ 0 :pcolumn-l ] ; 
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\   This  program  computes  the  LOFARGRAM  based  on  windowed  periodograms   and  displa%  ys  the  resi 

|id=no  of  points  in  xform,  step  =  shift,  tlen=transform  length 

I  [P] =lo farm (name, id, step, tlen) 

function  [P] =lofarm(name,  id, step, tlen) 

clear  pow 

%name=input ('    name   of    input    file    is  '); 

[drow, dcolumn] =size (name) ; 

if  drow~=l 

name=name . '  ; 
end 

len=length (name) ; 

%id=input('  No.  of  non-zero  data  points  in  transform       '  ) ; 
%step=mput  ('  shift  (step  size)  in  data  points      '); 
%tlen=input (' transform  length        '); 
i=fix ( (len-id) /step) +1; 
|win=hamming  (id)  .'  ; 
for  ic=l:i 

ppow=abs (f ft (name  (1+  (ic-1) *step: id+ (ic-1) *step) . *win,  tlen) ) ; 
pow (ic, : ) =ppow (1 : tlen/2) ; 
end 
P=pow; 
clear  name  drow  dcolumn  len  id  step  tlen  i  ic  win  ppow 
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9.  Appendix  A 

This  appendix  provides  3  figures  which  illustrate  the  effects  of 
level  selection  when  using  a  simple  (no-gray  scale  ,  no  color  ) 
printer.  All  three  figures  use  a  shift  of  4096  (i.e.  no  overlap  of 
data  between  successive  spectral  lines).  Fig  A.l,  A. 2  and  A. 3  show 
IPS,  CUM,  and  LOFAR  respectively.  The  number  of  quantization 
levels  is  different  on  each  figure.  The  label  on  each  figure  is 
self  explanatory. 
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